import json
import numpy as np
import pandas as pd
import scipy.optimize
import scipy.stats as st
import scipy.special
import numba
import tqdm
import biocircuits
import warnings
import iqplot
import bokeh.io
bokeh.io.output_notebook()
import panel as pn
pn.extension()
warnings.filterwarnings('ignore')
def style(p, autohide=False):
p.title.text_font="Helvetica"
p.title.text_font_size="16px"
p.title.align="center"
p.xaxis.axis_label_text_font="Helvetica"
p.yaxis.axis_label_text_font="Helvetica"
p.xaxis.axis_label_text_font_size="13px"
p.yaxis.axis_label_text_font_size="13px"
p.background_fill_alpha = 0
if autohide: p.toolbar.autohide=True
return p

$\hspace{15em}\mathrm{Nondimensionalizing... }\\[0.5em]$ \begin{align} \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}\tilde{t}} &= \beta_0 \cfrac{\tilde{c}(\tilde{m} m_d /k)^2}{1+(\tilde{m}m_d/k)^2} t_d - \gamma_c \tilde{c} t_d - \alpha_c m_d \tilde{m} \tilde{c} t_d \\[0.5em] \cfrac{\mathrm{d}\tilde{m}}{\mathrm{d}\tilde{t}} &= I_0\cfrac{t_d}{m_d} + \beta_1 c_d \tilde{c} \cfrac{t_d}{m_d} - \alpha_0 \cfrac{c_d \tilde{c} (m_d \tilde{m}/k)^2}{1+(m_d \tilde{m}/k)^2}\cfrac{t_d}{m_d} - \gamma_m \tilde{m} t_d \\[1.5em] \end{align}
$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{m_d = k, \hspace{0.5em} t_d = 1/\gamma_m, \hspace{0.5em} c_d = \cfrac{\gamma_m k}{\beta_1}}\\[0.1em] \end{align} \begin{align} \\[0.5em] \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}\tilde{t}} &= \cfrac{\beta_0}{\gamma_m} \cfrac{\tilde{c}\tilde{m}^2}{1+\tilde{m}^2} - \cfrac{\gamma_c}{\gamma_m} \tilde{c} - \cfrac{\alpha_c k}{\gamma_m}\tilde{m}\tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{m}}{\mathrm{d}\tilde{t}} &= \cfrac{I_0}{k \gamma_m} + \cfrac{\beta_1}{\gamma_m k} c_d \tilde{c} - \cfrac{a_0}{\gamma_m k}c_d \tilde{c} \cfrac{\tilde{m}^2}{1+\tilde{m}^2} - \tilde{m}\\[1.5em] \end{align}
$\hspace{15em}\mathrm{Setting... }$ \begin{align} \boxed{\tilde{\beta_0} = \beta_0/\gamma_m, \hspace{0.5em} \gamma = \gamma_c/\gamma_m, \hspace{0.5em} \tilde{\alpha_c} = \cfrac{\alpha_c k}{\gamma_m} \\[0.5em] \hspace{2em} \tilde{I_0} = \cfrac{I_0}{k \gamma_m}, \hspace{0.5em} \tilde{\alpha_0} = \alpha_0/\beta_1 } \end{align} \begin{align} \\[0.1em] \cfrac{\mathrm{d}\tilde{c}}{\mathrm{d}\tilde{t}} &= \tilde{\beta_0} \cfrac{\tilde{m}^2 \tilde{c}}{1+\tilde{m}^2} - \gamma \tilde{c} - \tilde{\alpha_c} \tilde{m} \tilde{c} \\[0.5em] \cfrac{\mathrm{d}\tilde{m}}{\mathrm{d}\tilde{t}} &= \tilde{I_0} + \tilde{c} - \tilde{\alpha_0} \tilde{c} \cfrac{\tilde{m}^2}{1+\tilde{m}^2} - \tilde{m}\\[1.5em] \end{align} $\hspace{15em}\mathrm{Rewriting... }$ \begin{align} \boxed{ \cfrac{\mathrm{d}c}{\mathrm{d}t} = {\beta_0} \cfrac{m^2 c}{1+m^2} - \gamma c - \alpha_c m c \\[0.5em] \cfrac{\mathrm{d}m}{\mathrm{d}t} = I_0 + c - \alpha_0 c \cfrac{m^2}{1+m^2} - m } \end{align}
No matter what the concentration of $m$, if $c = 0$, then c is not changing. Thus in the $c-m$ plane, this is a vertical line at $c = 0$, and $m$ takes on the entire range of real (positive) numbers. Note that when we plug in $c = 0$ to dm/dt, then we get $m = I_0$.
When $\Delta > 0$: \begin{align} \alpha_c^2 [(\gamma - \beta_0)^2 - 27 \gamma^2 18 \gamma (\gamma - \beta_0)] - 4 \gamma (\gamma - \beta_0)^3 - 4 \alpha_c^4 &> 0 \\[0.5em] \alpha_c^2 [(\gamma-\beta_0)(\gamma-\beta_0+18\gamma)-27\gamma^2] - 4 \gamma (\gamma - \beta_0)^3 - 4 \alpha_c^4 &> 0 \\[0.5em] \alpha_c^2 [-8\gamma^2 + \beta_0^2 - 20 \gamma\beta_0] - 4 \gamma (\gamma - \beta_0)^3 - 4 \alpha_c^4 &> 0 \\[0.5em] \boxed{\alpha_c^2 [\beta_0^2 - 20 \gamma\beta_0 - 8\gamma^2] + 4 \gamma (\beta_0 - \gamma)^3 - 4 \alpha_c^4 > 0}\\[0.5em] \end{align} I believe there is a typo in the problem statement, (the last term $4\alpha_c^4$ should have a power of 4 and not 2). I could be wrong though (I pulled the determinant equation off of Wikipedia), but I suspect it's a typo otherwise the squared term would be grouped together and simplified.
Rewriting \begin{align} \cfrac{\mathrm{d}m}{\mathrm{d}t} = 0 &= I_0 + c - \alpha_0 c \cfrac{m^2}{1+m^2} - m \\[0.5em] &0 = I_0(1+m^2) + c(1+m^2) - \alpha_0 c m^2 - m(1+m^2)\\[0.5em] &0 = -m^3 + I_0 m^2 - m + I_0 + c(1 + m^2 - \alpha_0 m^2)\\[0.5em] &c = \cfrac{m^3 - I_0 m^2 + m - I_0}{1+m^2(1-\alpha_0)}\\[0.5em] &\boxed{c = \cfrac{(m-I_0)(1+m^2)}{1+m^2(1-\alpha_0)}}\\[0.5em] \end{align}
This function is undefined when the denominator approaches zero, so: $$m_{st}^2 \neq \cfrac{1}{\alpha_0-1}$$ Furthermore, for c > 0 requires either that both numerator and denominator are positive, or both are negative. Note that when $\alpha_0 \leq 1$, the denominator is always positive. So we can break up our cases like so:
When $\alpha_0$ is less than 1, there is a sharp lower bound for steady state $m$. $$ \boxed{m_{st} \geq I_o }$$
When $\alpha_0$is greater than 1, either the first inequality holds or the second. The steady state concentration of the cytokine is sandwiched between these two parameters. \begin{align} \boxed{ \hspace{3em} I_0 < m_{st} < \sqrt{\cfrac{1}{\alpha_0-1}} \\[0.5em] \sqrt{\cfrac{1}{\alpha_0-1}} < m_{st} < I_0 } \end{align}
Let's rewrite our original ODE and collect the coefficient of $c$ for some insight:
$\cfrac{\mathrm{d}m}{\mathrm{d}t} = I_o + c\left(1-\alpha_0\cfrac{m^2}{1+m^2}\right) - m $.
It now becomes more clear what the threshold at $\alpha_0 = 1$ means; $m$'s response to concentrations of $c$ is either capable of being a net positive or negative (the hill function activation term is bounded by one and reaches one for m sufficiently large). So when $\alpha_0 > 1$, dm/dt is cleared (taken in by cells) faster than the per cell rate at which it is being secreted (the original dimensionalized inequality looks at $\alpha_0$ vs. $\beta_1$).
Thus, if $\alpha_0$ is high, cytokine uptake is so high it sets what the upper bound of what $m_{st}$ could ever be when $I_0$, the basal level of production, is small. Conversely, when $I_0$ is relatively high, $m$ can recover from uptake by the cells and is relatively unaffected. But because the cells are also secreting cytokine, there is an inherent lower bound within $\alpha_0$.
Analytically, we can think of $m_{st} (\alpha_c, \gamma, \beta_0)$, and $c_{st}(m_{st}, I_0, \alpha_0)$, where both steady states are determined purely by the adjacent parameters. From a design standpoint, if we wanted two stable states we want to have an $I_o$ and $\alpha_0$ that satisfies the relation above.
beta_o = 9.78
gamma = 2.14
Io = 24.67
alpha_c = 1.06
alpha_o = 11.21
c0_st = 0
m0_st = Io
_roots = np.roots([alpha_c, gamma-beta_o, alpha_c, gamma])
m1_st, m2_st = [_root for _root in _roots if _root >= 0]
find_c_cubic = lambda m : (m**3 - Io*m**2 + m - Io) / (1 + (m**2)*(1-alpha_o))
c1_st, c2_st = find_c_cubic(m1_st), find_c_cubic(m2_st)
print(" (m, c)")
print("({:.4}, {:})".format(m0_st, c0_st))
print("({:.4}, {:.3})".format(m1_st, c1_st))
print("({:.3}, {:.4})".format(m2_st, c2_st))
(m, c) (24.67, 0) (7.024, 1.77) (0.636, 10.8)
color_dc_dt = "#9fc0c1"
color_dm_dt = "#ecbac5"
color_fixed = "#65042d"
params_c = (alpha_c, beta_o, gamma)
params_m = (alpha_o, Io)
args = (alpha_o, alpha_c, beta_o, gamma, Io)
fps = np.array([[c0_st, m0_st], [c1_st, m1_st], [c2_st, m2_st]])
LOG_toggle = pn.widgets.Toggle(name="LOG", width=300)
@pn.depends(LOG_toggle.param.value)
def plotter_log_nullcline(LOG):
x_axis_type = "log" if LOG else "linear"
y_axis_type = "log" if LOG else "linear"
x_axis_label = "log(c)" if LOG else "c"
y_axis_label = "log(m)" if LOG else "m"
x_range = (1e-1, 1e3) if LOG else (-1, 16)
y_range = (1e-2, 1e3) if LOG else (-1, 29)
c_range = (-5, 4) if LOG else (-1.5, 20)
m_range = (-7, 7) if LOG else (-0.28, 35)
p = bokeh.plotting.figure(
height=375,
width=550,
title="Nullclines on the 𝑐−𝑚 plane",
x_axis_label=x_axis_label,
y_axis_label=y_axis_label,
x_range=x_range,
y_range=y_range,
x_axis_type=x_axis_type,
y_axis_type=y_axis_type,
)
# .... NULLCLINES ....
lw = 3.5
if LOG:
c_space = np.logspace(c_range[0], c_range[1], 1000)
m_space = np.logspace(m_range[0], m_range[1], 1000)
else:
c_space = np.linspace(c_range[0], c_range[1], 1000)
m_space = np.linspace(m_range[0], m_range[1], 1000)
c_zero = 0 * m_space
c_cubic = find_c_cubic(m_space)
p.line(c_cubic, m_space, line_width=lw, line_color=color_dm_dt)
p.line(c_zero, m_space, line_width=lw, line_color=color_dc_dt)
p.line(c_space, m1_st, line_width=lw, line_color=color_dc_dt)
p.line(c_space, m2_st, line_width=lw, line_color=color_dc_dt)
# .... FIXED POINTS ....
# unstable unfilled hole
p.circle(
c1_st, m1_st,
size=12,
color="white",
line_color="black",
line_width=4
)
# pizazz wedges
spin = 0
start_angles = np.linspace(0, 2*np.pi, 12)
end_angles = start_angles + 0.1
in_radius = 0.8 if LOG else 0.5
out_radius = 0.3 if LOG else 0.2
for start_angle, end_angle in zip(start_angles, end_angles):
p.annular_wedge(
x=c1_st,
y=m1_st,
fill_color="black",
line_color="black",
inner_radius=in_radius,
outer_radius=in_radius+out_radius,
start_angle=start_angle+spin,
end_angle=end_angle+spin,
)
# stable filled hole
p.circle(c0_st, m0_st, size=16, color="black", line_color="white", line_width=1.8)
p.circle(c2_st, m2_st, size=16, color="black", line_color="white", line_width=1.8)
# .... LEGEND ....
legend = bokeh.models.Legend(items=[
("stable fp", [p.circle(color="black", size=15)]),
("unstable fp", [p.circle(color="white", size=15, line_color="black", line_width=2)]),
("dc/dt nullclines", [p.line(color=color_dc_dt, line_width=lw)]),
("dm/dt nullcline", [p.line(color=color_dm_dt, line_width=lw)]),
], location="center")
p.add_layout(legend, 'right')
return style(p, autohide=True)
pn.Column(plotter_log_nullcline, pn.Row(pn.Spacer(width=50), LOG_toggle))
The log toggle does not add much, it just helps zoom in on the bottom right fixed point. (note that $c = 0 \rightarrow \log c = -\infty$, so the logged plot only shows one stable fixed point and is limited in this sense)
Let's rewrite our original ODE and collect the coefficient of $c$ for some insight: $\cfrac{\mathrm{d}m}{\mathrm{d}t} = I_o + c\left(1-\alpha_0\cfrac{m^2}{1+m^2}\right) - m $.
It now becomes more clear what the threshold at $\alpha_0 = 1$ means; $m$'s response to concentrations of $c$ is either capable of being a net positive or negative (the hill function activation term is bounded by one and reaches one for m sufficiently large). So when $\alpha_0 > 1$, dm/dt is cleared (taken in by cells) faster than the per cell rate at which it is being secreted (the original dimensionalized inequality looks at $\alpha_0$ vs. $\beta_1$).
Thus, if $\alpha_0$ is high, cytokine uptake is so high it sets what the upper bound of what $m_{st}$ could ever be when $I_0$, the basal level of production, is small. Conversely, when $I_0$ is relatively high, $m$ can recover from uptake by the cells and is relatively unaffected. But because the cells are also secreting cytokine, there is an inherent lower bound within $\alpha_0$. From a design standpoint Analytically, we can think of $m_{st} (\alpha_c, \gamma, \beta_0)$, and $c_{st}(m_{st}, I_0, \alpha_0)$, where both steady states are determined purely by the adjacent parameters. From a design standpoint, if we wanted two stable states we want to have an $I_o$ and $\alpha_0$ that satisfies the relation above.
Linear stability matrix \begin{align} A = \begin{pmatrix} \beta_0 \cfrac{m^2}{1+m^2} - \gamma - \alpha_c m, & \beta_0 \cfrac{2cm}{(1+m^2)^2} - \alpha_c c\\ 1 - \alpha_0 \cfrac{m^2}{1+m^2}, & - \alpha_0 \cfrac{2cm}{(1+m^2)^2} - 1 \end{pmatrix} \end{align}
def dc_dt(c, m, alpha_c, beta_o, gamma):
term1 = beta_o * m**2*c / (1 + m**2)
term2 = - gamma * c
term3 = - alpha_c * m * c
return term1 + term2 + term3
def dm_dt(c, m, alpha_o, Io):
term1 = Io + c
term2 = - alpha_o * m**2 * c / (1 + m**2)
term3 = - m
return term1 + term2 + term3
def ode_rhs(x, t, alpha_o, alpha_c, beta_o, gamma, Io):
c, m = x
dm_dt = Io + c - alpha_o * c * m**2/(1+m**2) - m
dc_dt = beta_o * m**2*c/(1+m**2) - gamma*c - alpha_c * m * c
return np.array([dc_dt, dm_dt])
def lin_stab_matrix(c, m, alpha_o, alpha_c, beta_o, gamma, Io):
e00 = beta_o * m**2 / (1+m**2) - gamma - alpha_c * m # df(c)/dc
e11 = - alpha_o * c * 2*m/((1+m**2)**2) - 1 # df(m)/dm
e01 = beta_o*2*c*m/((1+m**2)**2) - alpha_c * c # df(c)/dm
e10 = 1 - alpha_o * m**2/(1+m**2) # df(m)/dc
A = np.array([[e00, e01], [e10, e11]])
return A
def lin_stab(c, m, alpha_o, alpha_c, beta_o, gamma, Io):
A = lin_stab_matrix(c, m, alpha_o, alpha_c, beta_o, gamma, Io)
return np.linalg.eig(A)
def plot_separatrix(
c,
m,
c_range,
m_range,
p,
args,
t_max=50,
eps=1e-6,
color="#d95f02",
line_width=2,
line_dash="solid",
log=False,
return_separatrix=False,
):
# Negative time function to integrate to compute separatrix
def rhs(x, t):
c, m = x
# Stop integrating if we get the edge of where we want to integrate
if c_range[0] < c < c_range[1] and m_range[0] < m < m_range[1]:
return -ode_rhs(x, t, *args)
else:
return np.array([0, 0])
evals, evecs = lin_stab(c, m, *args)
evec = evecs[:, evals < 0].flatten() # Get eigenvector corresponding to attraction
t = np.linspace(0, t_max, 1000)
# Build upper right branch of separatrix
x0 = np.array([c, m]) + eps * evec
x_upper = scipy.integrate.odeint(rhs, x0, t)
# Build lower left branch of separatrix
x0 = np.array([c, m]) - eps * evec
x_lower = scipy.integrate.odeint(rhs, x0, t)
# Concatenate, reversing lower so points are sequential
sep_c = np.concatenate((x_lower[::-1, 0], x_upper[:, 0]))
sep_m = np.concatenate((x_lower[::-1, 1], x_upper[:, 1]))
if return_separatrix: return sep_c, sep_m
# Plot
if log:
p.line(np.log(sep_c), np.log(sep_m), color=color, line_width=line_width, line_dash=line_dash)
else:
p.line(sep_c, sep_m, color=color, line_width=line_width, line_dash=line_dash)
return p
Since this code is a direct copy-paste of what's above, except the phase portrait which is really just calling biocircuits, I've abstracted the code away into a module called utilities.py. Let's call the function and build our plot! (Note that to view the code of the function in notebook can run the following line: utilities.part_ef_plotter??)
import utilities
c_unstable, m_unstable = fps[1]
c_range=(-1.5, 13)
m_range=(-0.28, 28)
sep_c, sep_m = plot_separatrix(
c_unstable,
m_unstable,
c_range,
m_range,
0,
args,
return_separatrix=True
)
p = utilities.part_ef_plotter(sep_c[930:1200], sep_m[930:1200])
bokeh.io.show(p)
color_M = "#db98a5"
color_C = "#9ab5b6"
q = bokeh.plotting.figure(
height=360,
width=600,
title="Above and Below the Separatrix",
x_axis_label="dimensionless time",
y_axis_label="[ ]"
)
t = np.linspace(0, 7, 1000)
CMo_below = (1.8, 6.9)
CMo_above = (1.8, 7.2)
for CMo, line_dash in zip([CMo_below, CMo_above], ["solid", "dotdash"]):
_CM = scipy.integrate.odeint(ode_rhs, CMo, t, args=args)
C, M = _CM.T
q.line(t, M, line_color=color_M, line_width=4, line_dash=line_dash)
q.line(t, C, line_color=color_C, line_width=4, line_dash=line_dash)
legend = bokeh.models.Legend(items=[
("m above sep", [q.line(line_color=color_M, line_width=3, line_dash="dotdash")]),
("c below sep", [q.line(line_color=color_C, line_width=3, line_dash="solid")]),
("m below sep", [q.line(line_color=color_M, line_width=3, line_dash="solid")]),
("c above sep", [q.line(line_color=color_C, line_width=3, line_dash="dotdash")]),
], location="center")
q.add_layout(legend, 'right')
bokeh.io.show(style(q))
So we see above the separatrix, we converge to the concentrations (c, m) of the stable fixed point (0, 24.7), and below the separatrix, we converge to the point (10.8, 0.64). Looking at the shape of the separatrix, we see that to reach the homeostatic stable state (so that cytosine concentration doesn't just collapse), the initial m generally has to be below ~2-3x the initial c.
The large divisive basins provide a threshold for what we colloquially refer to as homeostasis, where perturbations within the basin below the separatrix demonstrates a robustness to the initial starting concentrations, only sensitive along the separatrix itself (as demonstrated in the plot above).
Otherwise the system's steady state remains robust to a wide range of initial conditions (shown below). Here, we take our separatrix, then shift sep_m by 1 below the separatrix line, and plot all the traces of this trajectory. Since the code is really just the same as what's shown above, the function is stored in utilities.py as part_g_plotter().
utilities.part_g_plotter(sep_c, sep_m)
I was going to try a stochastic simulation next, but I got distracted and found this lil easter egg in the comments of the problem statement!
We are cutting this part of the problem just to keep the homework from getting too long.
h) Play with parameter values and investigate conditions for bistability. You already have much the machinery in place to do this, since you can already computed the fixed points. For bistability, we must have three total fixed points (two in addition to the $C = 0$ fixed point that always exists). You can use this as a criterion to find ranges of parameter values for which bistability exists. With respect to which parameter values is bistability most robust? How about least robust?
I abstracted a lot of the code away into a module utilities.py. The major addition is defining a function draw_fp() that will responsively plot different glyphs for different stabilities, defining a dashboard-creator param_plotter() that will plot our nullclines as we toggle our parameters. Here is the result:
import utilities
_param_plotter = utilities.param_plotter()
_param_plotter.servable()
Note that as soon as roots become complex, the points disappear (try-except-pass to the rescue!)
MORE ANALYSIS
So we make the following conclusion from our mathematical analysis in parts (a) - (g) and the visualizations above. The existence of $m_{st}$ depends on the existence of positive real roots in the cubic that is parametrized by three parameters: $(\alpha_c, \gamma, \beta_0)$. We saw in part (b) that a necessary condition is that $\gamma < \beta_0$ and that the discriminant, a quadratic in $\alpha_c^2$, be positive. Once the root of $m_{st}$ is found, $c_{st}$ can be found by plugging in $m_{st}$ to an equation parametrized by $(\alpha_0, I_0)$. The positive nature of $c_{st}$ is then dependent on the inequalities we set up in part (d.iv).
The authors of the paper provided values that we used for parts (a)-(g). Now, we will broaden our search:
classifier() in utilities.py (please check that I did this correctly, as my entire analysis rests on this function working the way it should) that takes in five parameters, and returns 0 for no fixed point, 1 for an unstable fixed point, and 2 for a stable fixed point. .... Although this function isn't called directly, but through utilities.py, I've reproduced it below. (pseudocode: first check length of positive real $m_{st}$ array. If zero, return zero. then iterate through all fixed points. if one has both eigenvalues as negative, and positive c from the Io and alpha conditions, return 2. else, return 1). Note that I do not handle the case where the eigenvalue is zero..json files, titled abg.json for the $(\alpha_c, \beta_0, \gamma)$ exploration and aig.json for the $(I_0, \alpha_0, \gamma)$ exploration. The code is provided in Appendix B (at the end of this notebook). beta_o = 9.78
gamma = 2.14
Io = 24.67
alpha_c = 1.06
alpha_o = 11.21
def classifier(alpha_o, alpha_c, beta_o, gamma, Io):
"""
Returns 0: positive real m, positive c NOT found
Returns 1: positive real m, positive c, unstable fp
Returns 2: positive real m, positive c, stable fp
"""
ms_st = []
_roots = np.roots([alpha_c, gamma-beta_o, alpha_c, gamma])
for _root in _roots:
if (np.imag(_root) == 0.0) & (np.real(_root) > 0.0):
ms_st.append(_root)
cs_st = [find_c_cubic(m, alpha_o, Io) for m in ms_st]
if len(ms_st) == 0:
return 0
args = (alpha_o, alpha_c, beta_o, gamma, Io)
for m, c in zip(ms_st, cs_st):
evals, evecs = lin_stab(c, m, *args)
if (np.all(evals < 0)) and (c > 0):
return 2
return 1
## .... LOADING IN JSON FILES ....
with open('abg.json', 'r') as f:
d_abg = json.load(f)
with open('aig.json', 'r') as f:
d_aig = json.load(f)
Even with the json files, the slider is really slow on my computer (150x150 grid for 2 plots with a slider that takes on 20 values means 900,000 points in total) so I tried turning it into a movie—hopefully it works!
from IPython.display import Video
Video("all_out.mp4", width=900)
Here is the dashboard:
abig_dashboard = utilities.stab_abig_plotter(d_abg, d_aig)
abig_dashboard.servable()
So we see that we do not want $\gamma$ to be too high, otherwise we won't be able to find an $I_0$ or $\alpha_0$ at a given $\alpha_c$, $\beta_0$. But once we satisfy that condition, we can have a stable fixed point when both $I_0$ and $\alpha_0$ are of similar orders of magnitude. The block-like shape emerges given the chained inequalities we set up earlier, where our root is satisfied either with $I_0$ as a lower bound (upper right block), or as an upper bound (lower left block), and vice veras for $\alpha_0$.
To answer the original question posed in the prompt, the system appears most robust to $Io, \alpha_0$, given the broad ranges they can take on. In regions of small $\alpha_c$, $\beta_0$ can take a range of values (big red block at $\gamma = 5.5$). I am hesitant to say that the system is most sensitive to $\gamma$, since that is the parameter I have chosen as my slider; it could easily be the case we lose stability in the right plot when looking along other dimensions. However, it is certainly the relation between $\beta_0$ and $\gamma$ that results in greater sensitvity, for we see that $alpha_c$ can take on values from $10^{-3} - 10^2$, whereas $\beta$ and $\gamma$ are both sandwiched for the most part between $10^1 - 10^2$. (note that although the left plot is a square, the scale of the axes are not, so the y-chunk is narrower than the x-chunk in actual parameter space).
Thank you much for reading!! This is all for 7.1.
Note: Is it ok that I am putting functions in a separate module? I just find the scrolling to be so overwhelming that I want the focus to be on my analysis and the visuals instead of the enormous amount of semi-repetitive code. Please let me know if this is counterproductive, as I can understand why it is important to be able to look at and interact with the code to see if it is functioning as it claims. That being said, I've included an "Appendix" for my functions in utilities, and the code that was used to create the .json files
## .... UNCOMMENT TO LOOK AT FUNCTIONS ....
## ----------------------------------------
# utilities.draw_fp??
# utilities.param_plotter??
# utilities.classifier??
Here is the code that was used to create the abg.json and aig.json files.
# ****** NOTE: a: alpha_c, b: beta_o ******
# ng = 20
# states = np.empty((ng, na, nb))
# gamma_range = np.logspace(-1, 2, ng)
# for k, gamma in enumerate(tqdm(gamma_range)):
# na, nb = 150, 150
# a = np.logspace(-5, 5, na)
# b = np.logspace(1, 4, nb)
# aa, bb = np.meshgrid(a, b, indexing='xy')
# for i in range(na):
# for j in range(nb):
# _a, _b = aa[i, j], bb[i, j]
# states[k, i, j] = classifier(alpha_o, _a, _b, gamma, Io)
# color_dict = {0: '#eba8b5', 1:'#9fc0c1', 2:'#65042d'}
# _df_alphas, _df_betas, _df_colors, _df_gammas = [], [], [], []
# for k, gamma in enumerate(gamma_range):
# for i in range(na):
# for j in range(nb):
# alpha, beta = aa[i, j], bb[i, j]
# color = color_dict[states[k, i, j]]
# _df_alphas.append(alpha)
# _df_betas.append(beta)
# _df_gammas.append(gamma)
# _df_colors.append(color)
# _df_alphas = [float(_) for _ in _df_alphas]
# _df_betas = [float(_) for _ in _df_betas]
# _df_gammas = [float(_) for _ in _df_gammas]
# d_abg = {'alpha':_df_alphas, 'beta':_df_betas, 'gamma':_df_gammas, 'color':_df_colors}
# df = pd.DataFrame(d_abg)
# import json
# ## .... SAVING AS JSON FILE ....
# with open('abg.json', 'w') as f:
# json.dump(d_abg, f)
## ****** NOTE: a: alpah_o, b: Io ******
# ng = 20
# states = np.empty((ng, na, nb))
# gamma_range = np.logspace(-1, 2, ng)
# for k, gamma in enumerate(tqdm(gamma_range)):
# na, nb = 150, 150
# a = np.logspace(-2, 2, na)
# b = np.logspace(-5, 5, nb)
# aa, bb = np.meshgrid(a, b, indexing='xy')
# for i in range(na):
# for j in range(nb):
# _a, _b = aa[i, j], bb[i, j]
# states[k, i, j] = classifier(_a, alpha_c, beta_o, gamma, _b)
# _df_alphas, _df_betas, _df_colors, _df_gammas = [], [], [], []
# for k, gamma in enumerate(gamma_range):
# for i in range(na):
# for j in range(nb):
# alpha, beta = aa[i, j], bb[i, j]
# color = color_dict[states[k, i, j]]
# _df_alphas.append(alpha)
# _df_betas.append(beta)
# _df_gammas.append(gamma)
# _df_colors.append(color)
# _df_alphas = [float(_) for _ in _df_alphas]
# _df_betas = [float(_) for _ in _df_betas]
# _df_gammas = [float(_) for _ in _df_gammas]
# d_aig = {'alpha_o':_df_alphas, 'Io':_df_betas, 'gamma':_df_gammas, 'color':_df_colors}
# df = pd.DataFrame(d_aig)
# ## .... SAVING AS JSON FILE ....
# with open('aig.json', 'w') as f:
# json.dump(d_aig, f)
~ Poem? ~
Rise and shine, my Cytokine.
Wine and dine,
down my spine,
align your signs—
for bistability is divine!